function foreground_detection(path, filename, frames, initframes)
% FOREGROUND_DETECTION is the main routine based on [1].
% 
% Project:  Foreground detection [Bildfolgen LU]@CVL/TUWIEN
% Author:   Gruppe 10
% Date:     03/2011

% [1] Elhefnawy, W.R., Selim, G., Ghoniemy, S.: Robust Hybrid Foreground Detection 
%       and Adaptive Background Model. International Conference on Information 
%       Science and Applications (ICISA), pp. 1-8 (2010) 

% clc;
% close all;
% clear all;
% path = 'S3-A320-2\\';
% filename = 'S3A3202'; %, 'S21-2', 'S36-A319-8'};
% frames = 50;
% initframes = 50; 



%ALPHA_FG = .001;    
ALPHA_BG = .05;%.15;
FILE_EXT = 'jpeg';
% h = fspecial('gaussian', [5 5], .5);
se = strel('square', 10);

fprintf(1,'Gruppe 10 - FOREGROUND DETECTION:\n---------------------------------\n\n');
fprintf(1, 'Path: \t\t%s\nFilename: \t%s\nFrames: \t%i\nInitframes: %i\n\n',...
    path(1:end-1), filename, frames, initframes);

try 
    sz = size(imread(sprintf('%s%s_0000.%s', path, filename, FILE_EXT)));
    %sz([1 2]) = sz([1 2])/2;
catch exception
    fprintf(1, 'ERROR: File %s%s_0000.%s not found! (%s)', path, filename, FILE_EXT,...
        exception.identifier);
end

% initizialisation
fprintf(1,'1. Init background model:\n');
imgs = zeros([sz initframes]);
mu2 = zeros(sz);
var2 = zeros(sz);
var3 = zeros(sz);

fprintf(1, ' Determine background parameters...\n\n\n\n\n\n\n'); tic;
% for idx = initframes-1:-1:0
for idx = 0:initframes-1 
    imgname = sprintf('%s%s_%04d.%s', path, filename, idx, FILE_EXT);
%     imgs(:,:,:,idx+1) = imread(imgname); %imresize(imread(imgname), sz([1 2]));
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     mu2 = mu2 + double(imread(imgname));
    I = double(imread(imgname));
    xn_1 = I-mu2;
    mu2 = mu2 + xn_1/(idx+1);
    var2 = (idx*var2 + (I-mu2).*xn_1)/(idx+1);
    var3 = var3+xn_1.*(I-mu2);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %    imgs(:,:,:,idx+1) = rgb2i1i2i3(imgs(:,:,:,idx+1));
%     imgs(:,:,:,idx+1) = imfilter(imgs(:,:,:,idx+1), h, 'symmetric', 'conv');
%     fprintf(1, '\b\b\b\b\b\b\b%6.2f%%',(initframes-idx)/initframes*100);
    if(idx)
        fprintf(1, '%c', repmat(char(8), [1 length(imgname(numel(path)+1:end))+3]));
    end    
    fprintf(1, '\b\b\b\b\b\b\b%6.2f%% (%s)', (idx+1)/initframes*100, imgname(numel(path)+1:end));
end
fprintf(1, '%c', repmat(char(8), [1 length(imgname(numel(path)+1:end))+3]));
fprintf(1, '\t\t(%5.2f sec)\n', toc);
var3 = var3/initframes;


% fprintf(1, '  Determine background parameters...'); tic;
% [mu, vari] = initBGparam(imgs);
% fprintf(1, 'OK!\t(%5.2f sec)\n', toc);

mu = mu2;
vari = var3;
mu_dist = zeros(sz(1:2));
var_dist = zeros(sz(1:2));

fprintf(1, '\n2. Init shadow model:\n');
fprintf(1, '  Train shadow model...\n\n\n\n\n\n\n'); tic;
for idx = 0:initframes-1
    imgname = sprintf('%s%s_%04d.%s', path, filename, idx, FILE_EXT);
    tmp = double(imread(imgname))./mu;
    dist = arrayfun(@getDist, tmp(:,:,1), tmp(:,:,2), tmp(:,:,3));
%     delta = dist - mu_dist;
%     mu_dist = mu_dist + delta/(idx+1);
    var_dist = var_dist + dist.*dist;%delta.*dist; %(dist-mu_dist);
    fprintf(1, '\b\b\b\b\b\b\b%6.2f%%',(idx+1)/initframes*100);
end
clear imgs;
var_dist = var_dist./initframes;
var_dist = sqrt(mean(mean(var_dist)));
% likehood = arrayfun(@getLikelihood, mu_dist, var_dist); 
fprintf(1, '\t\t\t\t\t(%5.2f sec)\n', toc);
 


fprintf(1,'\n3. Foreground detection:\n');
fprintf(1, '  Compute foreground objects...\n\n\n\n\n\n\n'); tic;
for idx = 0:frames-1
    imgname = sprintf('%s%s_%04d.%s', path, filename, idx+initframes, FILE_EXT);
    %imgs(:,:,:,idx+1) = im2double(imread(imgname)); %imresize(imread(imgname), sz([1 2]));
    img = imread(imgname);
%     img = rgb2i1i2i3(img1);
%     img = imfilter(img, h, 'symmetric', 'conv');
    % foreground mask
    [fm dif] = getFGmask(img, mu, vari);    
    
    
    

    subplot(321)
    imshow(dif(:,:,1), []);
    title('R');
    subplot(323)
    imshow(dif(:,:,2), []);
    title('G');
    subplot(325)
    imshow(dif(:,:,3), []);
    title('B');
    fm = bwareaopen(fm, 40, 4);
    ifm = repmat(fm, [1 1 3]);
    rp = regionprops(fm, 'Area', 'BoundingBox', 'Centroid', 'MajorAxisLength',...
        'MinorAxisLength', 'Orientation');
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
    % SHADOW 
    sidx = find(ifm);
    delta = double(img(sidx))./mu(sidx);
    delta = reshape(delta, sum(sum(fm)), 3);
    dist = arrayfun(@getDist, delta(:,1), delta(:,2), delta(:,3));
    sm = normpdf(dist, 0, var_dist); %sqrt(var_dist(fm)))> 0.05;  %shadow mask
    sm = sm>quantile(sm, .50);
    fmidx = find(fm);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    %imshow(fm);
    out = img;
    %out(fm) = 255;
    %imshow(out)
    %hold on;
    edgm = zeros(size(fm));
    cedgm = zeros(size(fm));
    if(idx)
        fd = abs(double(img)-double(prev));
%         edgm = colfilt(fd(:,:,1), [3 3], 'sliding', @comp4Neigh) ...    % edge mask
%             | colfilt(fd(:,:,2), [3 3], 'sliding', @comp4Neigh) ...
%             | colfilt(fd(:,:,3), [3 3], 'sliding', @comp4Neigh);

       edgm = comp4Neigh(fd(:,:,1)) | comp4Neigh(fd(:,:,2)) | comp4Neigh(fd(:,:,3));   % edge mask
       ledgm = bwconncomp(bwareaopen(edgm, 40, 4), 8);   % labeled edge mask
       
       for i = 1:ledgm.NumObjects    
            cot = zeros(size(fm));  % core template object
            cot(ledgm.PixelIdxList{i}) = 1;          
            cedgm = cedgm | fillArea(cot);
       end
    end
    
    cbem = zeros(size(fm));   
    for bl = 1:numel(rp)    % all blobs
%         ellipse(rp(bl).MajorAxisLength/2, rp(bl).MinorAxisLength/2, rp(bl).Orientation*pi/180,...
%             rp(bl).Centroid(1), rp(bl).Centroid(2));
  %      rectangle('Position', rp(bl).BoundingBox, 'EdgeColor', 'g');
%         plot(rp(bl).Centroid(1), rp(bl).Centroid(2), 'cx'); % regionprops centroid
        
%         [y x] = find(lfm==bl);
%         mx = mean(x);
%         my = mean(y);
%         plot(mx,my, 'ko');    % own centroid(1st moment);
        
        el = getEllipse(rp(bl).MajorAxisLength/2, rp(bl).MinorAxisLength/2,... 
            rp(bl).Centroid(1), rp(bl).Centroid(2), rp(bl).Orientation*-1*pi/180); % own ellipse
        
        el(el(:,1)<1,1) = 1; el(el(:,1)>size(fm,2),1) = size(fm,2); 
        el(el(:,2)<1,2) = 1; el(el(:,2)>size(fm,1),2) = size(fm,1);
        bem = zeros(size(fm));   % bounding ellipse mask;
        bem(sub2ind(size(fm), round(el(:,2)), round(el(:,1)))) = 1;
        bem = fillArea(bem);     % filled bounding ellipse
  %      plot(el(:,1), el(:,2), 'c');
        cbem = cbem | bem;      % complete bounding ellipses
        %image(uint8(cbem)*255, 'AlphaData', double(cbem));
    end
    subplot(322)
    imshow(cedgm);
    subplot(324)
    out2 = out;
    %out2(fm&cedgm) = 255;
    out2(fmidx(sm)) = 255;
    imshow(out2);
    out(fm&cbem) = 255;
    subplot(326);
    imshow(out)
    %hold on;
    
    pause(.1);

%     % update FG
%     mu(ifm) = (1-ALPHA_FG)*mu(ifm) + ALPHA_FG*img(ifm);
%     sigma(ifm) = (1-ALPHA_FG)*sigma(ifm) + ALPHA_FG*dif(ifm);
    % update BG
    mu(~ifm) = (1-ALPHA_BG)*mu(~ifm) + ALPHA_BG*double(img(~ifm));
    vari(~ifm) = (1-ALPHA_BG)*vari(~ifm)+ ALPHA_BG*dif(~ifm).^2;
    
    
    
    if(idx)
        fprintf(1, '%c', repmat(char(8), [1 length(imgname(numel(path)+1:end))+3]));
    end    
    fprintf(1, '\b\b\b\b\b\b\b%6.2f%% (%s)', (idx+1)/frames*100, imgname(numel(path)+1:end));
    
    prev = img;     % previous frame 
    
end
fprintf(1, '%c', repmat(char(8), [1 length(imgname(numel(path)+1:end))+3]));
fprintf(1, '\t(%5.2f sec)\n', toc);


input('\nPress Enter to continue...\n');
%end


